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We investigate the occurrence of waterlike thermodynamic and dynamic anomalous behavior in a 
one dimensional lattice gas model. The system thermodynamics is obtained using the transfer matrix 
technique and anomalies on density and thermodynamic response functions are found. When the 
hydrogen bond (molecules separated by holes) is more attractive than the van der Waals interaction 
(molecules in contact) a transition between two fluid structures is found at null temperature and high 
pressure. This transition is analogous to a 'critical point' and intimately connects the anomalies in 
density and in thermodynamic response functions. Monte Carlo simulations were performed in the 
neighbourhood of this transition and used to calculate the self diffusion constant, which increases 
with density as in liquid water. 
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I. INTRODUCTION 

Liquid water has many properties which are recognized 
as being anomalous when compared to other nonbonded 
liquids with the same molecular size. As an example, 
the isothermal compressibility and the constant pressure 
heat capacity present minimum as a function of tempera- 
ture at ambient pressure. In addition, the thermal expan- 
sion coemcient is negative below 4°C, indicatingthat den- 
sity increases anomalously with temperature [H Q ■ Be- 
sides thermodynamic anomalies, water present dynamic 
anomalies. For temperatures below 283°C there is a re- 
gion where the self-diffusion constant increases as a func- 
tion of pressure [l|, H| . 

The thermodynamic anomalies seems to be inter- 
related and different thermodynamic scenarios were pro- 
posed to describe these relations Here we will focus 
our discussion on the second critical point scenario 
that is supported by the model studied here (as will be 
shown latter). According to this scenario, the diverging 
thermodynamic behavior of the response functions is as- 
sociated with the critical point arising from metastable 
liquid-liquid phase transition that occurs in the super- 
cooled regime, in an experimentally unaccessible temper- 
ature below the homogeneous nucleation temperature. 
This scenario was originally observed and proposed by 
means of computer simulations of atomistic models for 
liquid water [4[], but there are indirect experimental evi- 
dences for the existence of a liquid-liquid phase transition 
in supercooled water on amorphous [f| Q and confined 
water [7|. 

The relation between thermodynamic and dynamic 
anomalies have been subject of discussion in the liter- 
ature Using atomistic models of water, Errington 
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and Debenedetti found that the anomalously diffusive re- 
gion surrounds the region of density anomaly in the tem- 
perature vs. density plane Q. In addition, measures of 
translational and rotational structure also reveals unex- 
pected behavior, that could be connected to the anoma- 
lies on density and diffusion t 8|. These relations, usually 
called hierarchy of anomalies, are not restricted to wa- 
ter models and where observed in computer simulations 
of SiO ? \Hj [ill. BeF 2 [Hj], and core-softened models of 
fluids gflllHl||. 

With the aim of investigating the relation between 
thermodynamic and dynamic anomalies of water we pro- 
pose a one dimensional lattice model of water. Lattice 
and off-lattice one dimensional models of water were pro- 
posed by Ben- Nairn (T5T - [l7| . Bell[3], and others fToL T20jj . 
to investigate the so-called thermodynamic anomalies 
and even the unusual solvation behavior presented by 
water. Nevertheless, and as far as we know, the behav- 
ior of the equilibrium diffusion constant of these models 
were not investigated. 

The thermodynamics of our model is obtained exactly 
using transfer matrix technique and its dynamics is in- 
vestigated through Monte Carlo simulations. A watcrlikc 
anomalous behavior is found both on density and self dif- 
fusion constant, as observed in water [l| and atomistic 
models of water [1, [2l[ . 

This paper is organized as follows. In section |TT] we 
build the model's Hamiltonian and investigate its ground 
state. In section IIIII the thermodynamics of the model 
is analysed and discussed. Monte Carlo simulations are 
used to calculate the self-diffusion constant on section HVl 
and the relation between thermodynamics and kynetics 
is explored on section TlV Al Final remarks are made on 
the last section. 
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II. THE MODEL 

The model consists of a linear lattice whose sites can be 
either occupied by molecules or empty. Two interactions 
arc defined: a short-range van der Waals attraction be- 
tween nearest neighbours and a hydrogen bond between 
second nearest neighbours separated by holes. An occu- 
pation variable rjk is assigned to each site k to indicate the 
presence {rj k — 1) or the absence {rj k = 0) of a molecule. 
With these definitions the effective Hamiltonian in the 
grand canonical ensemble is written as: 

H = - [evdwViVi+i + £hbr)i(l - ^+1)^+2] 

i 

~P ^Vi, (!) 

i 

where e v d w > defines the strength of the van der Waals 
attraction, e^b > the strength of the hydrogen bond, 
and p is the chemical potential. 

In order to obtain a preliminary insight about the 
model we investigate the ground state by looking at three 
states putative states at null temperature: a gas (G), a 
bonded fluid (BF) and a dense fluid (DF) . The gas phase 
is an empty lattice which 'coexists' with the fluid phase 
of lower free energy at null pressure. The bonded fluid is 
a half filled and fully bonded lattice, while the dense fluid 
is a filled lattice without bonds. For a lattice with size 
L and periodic boundary conditions, the grand canonical 
free energies of the fluid phases are $bf = — (thb + p)L/2 
and $df = — {evdw + p)L. We will only consider pa- 
rameters that correspond to the hydrogen bond being 
more attractive than the van der Waals interaction, i.e., 
£/ifc < £vdw> m consistence with real water. In fact, this 
condition also ensures the predominance of a bonded 
fluid at low temperatures and pressures in the range 
< P < P c , with P c = ehb — tvdw At null tempera- 
ture, a transition between the BF and the DF happens 
at P c , as observed by Sadr-Lahijany et al. in a similar 
continous model [20| . Here we follow ref. [20] on calling 
this null temperature transition as the 'second critical 
point', in addition to the 'critical point' found at T = 
and P = 0. Further discussions about the 'second critical 
point' of this system will be made on section UlI Al 



III. THERMODYNAMICS 

We use the transfer matrix technique to exactly calcu- 
late the grand canonical partition function (more details 
can be found on appendix [X| of Hamiltonian ([1} as 




FIG. 1. Temperature of maximum density (TMD) in the p = 
Pl/tvdw vs. t = kT/e vdw phase diagram, for e hb /e vdw = 3/2. 
Note that the TMD starts in the low pressure 'critical point' 
at (0,0) and ends in the high temperature 'second critical 
point' at (0,0.5). 



mann constant. From ^(T,L,p) the density p becomes: 



P = 



where a = e^ e " d 



(A-l)[(aA + 6)(A-l) + l] 
(A-I) 2 (aA + 2&) + 3A-2 



(3) 



and b - 



It is shown on appendix 



A that the density can be written in terms of these quan- 
tities as: 

In addition to the density p, the densities (per site) 
of hydrogen bonds (phb) an d nearest neighbours (p nn ) 
are also important to relate the liquid structure to the 
anomalous behavior and are calculated on appendix A. 
The mean number of hydrogen bonds (nhb) and near- 
est neighbours (n nn ) per particle can be calculated from 
those expressions as n%b = Phb/ P and n nn — p nn /p. 

The thermodynamic response functions are important 
to characterize waterlike behavior, particularly because 
they indicate the possibility of a critical behavior on the 
supercooled regime (as is expected to occur in liquid wa- 
ter). Here we investigate the thermal expansion coeffi- 
cient (a) and isothermal compressibility (k), defined as: 



1 ( dp 

ct = t; — 

P \dT 

R= l(^P_ 
p\dP /7 



(4a) 
(4b) 



Expressions for these functions are lengthy and they will 
not be presented here. 



~(T,L,p) = \ I 



(2) 



where T, L, p and P are, respectively, the thermody- 
namic variables temperature, system's length, chemical 
potential and pressure, and j3 = 1/kT, with k the Boltz- 



A. Low temperature limit 

Let us note that for P = expression results in a 
null density p = 0. Nevertheless, at null temperature the 



3 



TABLE I. Value of various quantities while approaching the 
transition at (T c = 0, P c = e^b — tvdw) through different limits. 
Only the temperature dependence is presented for a and kr- 





P^P- 


P = Pc 


P^P+ 


p 


1/2 


2/3 


1 


phb 


1/2 


1/3 





Pnn 





1/2 


1 


a ~ 


-1/T 





1/T 


kx ~ 


1/T 


1/T 


1/T 



limiting value of the density depends on how the point 
(T = 0,P = 0) is approached, due the 'transition' oc- 
curing at this point. While approached from above, the 
density is compatible with the BF ground state (see the 
second column of Tab. |IJ. In addition, both a and Ut di- 
verges while approaching this 'critical point' (not shown). 
An equivalent behavior was previously obtained by Sadr- 
Lahijany et al. in a similar continous model [20]. Even 
more interesting is the low-temperature behaviour of the 
system in the neighbourhood of the second critical point 
(T = 0,P = P c ). In what follows we analyse this be- 
havior focusing on density p and its derivative a. From 
now on we assume P > 0. Using this one obtains the low 
temperature limit of eq. (j3|) as 



aX + b 
aA + 26' 



(5) 



where it was assumed that A ^> 1 at small temperatures. 
By applying this expression to the point (T = 0, P = 0) 
it follows that p — ► 1/2 while the point is approached 
from above, as discussed before. 



Now, we note that A r 



b/a in the neighbour- 



hood of the fluid-fluid transition at 



T — T c — 0, P — P c — Chb ~ tvdw 



and define 



P = Pc + A, 



A = KS, 5 = e /3A , 



in order to simplify eq. ([5]) . While approaching this tran- 
sition it can be found that 8 — > for A > 0, 5 = 1 for 
A = and <5 _1 — > when A < 0. Using these definitions, 
eq. ([5]) is written as 



1 + 8 

2 + 8' 



(6) 



The low-temperature limits of the density and the deriva- 
tives of © are obtained in a straightforward manner 
(see Tab. UJ). It is particularly interesting to note that 
the values expected from section |TT] are reobtained when 
the 'critical' pressure is approached from below or from 
above, see table HI Nevertheless, when null temperature 
limit is approached with P = P c one finds p — > 2/3, which 
is different from both BF and DF. This is an indication of 



a new liquid structure, which can be represented through 
a serie of two neighbouring water molecules surrounded 
by empty sites. By comparing the expected grand canon- 
ical free energies, it was verified a posteriori that this 
state coexists with both BF and DF exactly at P c . The 
existence of such a hidden state make the analysis of the 
ground state limit more complex. Nevertheless, a sim- 
ple derivation similar to the one presented above, was 
used for obtaining the low-temperature limiting values of 
the densities of hydrogen bonds (|A9I) and nearest neigh- 
bours (|A10|) presented on table fl] 

We finish this analysis noting that the behavior of the 
thermal expansion coefficient in the neighbourhood of the 
fluid-fluid phase transition changes qualitatively depend- 
ing on how the 'second critical point' is approached. For 
a linear approach A = kTtanO we find a cx tan0/T, 
which is reasonable since the density anomalous increase 
with temperature in region dominated by BF fluid struc- 
ture, for A < or equivalently P < P c . 



B. Liquid structure 

In what follows we analyse the thermodynamic liquid 
structure in the pressure vs. temperature phase diagram 
for the parameter ehb/^vdw = 3/2. Similar features are 
found for other parameter values, provided that e vc i w < 

Fig. [1] shows that the line of temperature of maximum 
density (TMD) starts at T = exactly at the BF-DF 
transition, which has the peculiar property that 



a(0,P c ) 



lim a(T, P) -> 0. 

T->0,P=P c 



This line reaches its maximum temperature at t = 
kT/e v dw ~ 0.388(1) and then its temperature decreases 
with pressure until reaching the G-BF transition located 
at (T c — 0,P C = 0). As discussed previously by Sadr- 
Lahijany et al. [20], the G-BF transition is a remanent 
of a gas-liquid phase transition of a system with higher 
dimensionality, and can be found in simple ID fluids 
with attractive short-range interactions. The very fact 
that the TMD connects two transitions (G-BF and BF- 
DF) seems to be a peculiar property of some ID mod- 
els, since a similar feature also occurs in the ID lattice 
model proposed by Bell [l8j and for some parameters in 
reference [201 ] . but not int the continuous models studied 

in [H & 

In Fig. [2] the behavior of density, number of hydro- 
gen bonds and number of nearest neighbours per particle 
is compared as a function of temperature and pressure. 
Figs. [21(a), (c) and (e) show that the anomalous increase 
in density with temperature (for P < P c ) is followed by 
a decrease on the number of hydrogen bonds and by an 
increase on the number of first nearest neighbours. For 
P = P c , not only a ~ is persistent for a wide range 
of temperatures, but (dnhb/dT)p w (dn nn /dT)p w is 
also found in the same region. Figs. [2£b), (d) and (f) 
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FIG. 2. The temperature dependence of p (a), nub an d n nn is shown on the left side for p 
0.55 (solid). The pressure dependence of the same quantities is shown on right side for t 
(dashdot). Units are the same used in Fig. [T] 



= 0.45 (dashdot), 0.50 (dash) and 
= 0.05 (solid), 0.1 (dash) and 0.3 
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FIG. 3. Isothermal compressibility function of temper- 

ature at same pressures as in Fig. [2] p — 0.45 (solid), 0.50 
(dashed) and 0.55 (dashdot). 



FIG. 4. Thermal expansion coefficient a as a function of 
pressure at temperatures t = 0.02 (solid), 0.05 (dashed) and 
0.1 (dashdot). 



show that the continuos transition between BF and DF 
progressively becomes more abrupt as the temperature is 
lowered. From these figures, it is also evident that p, nhb 
and n nn are essentialy not changing with temperature 
at P = P c , since different isotherms cross at this point. 
These results indicate that the intermediate fluid struc- 
ture discussed in the previous section is well populated in 
a borderline separating a BF-like region from a DF-like 
region. It is possible that this population occurs because 
the intermediate structure is an unavoidable intermediate 
step for changing between the BF and DF liquid struc- 
tures. A last comment about Figs. [2] (b), (d) and (f) is 
approapriate at this point: it is evident from the high 
temperature curves that a fast transition between the G 
and BF structures takes place at low pressures. This is 
associated with the low-T and low-P part of the TMD, 
which happens with a fast but small increase in density 
with temperature, as also observed in Figs. 5 and 6 from 
reference fl8j |. 

The isothermal compressibility k is presented in Fig. [3] 
as a function of temperature, at the same pressures as 
Figs. Eta), (c) and (e). Note that k increases with de- 
creasing temperature in the neighbourhood of of the BF- 
DF transition, and that it diverges when approaching 
this transition with P = P c . The diverging behavior of 
a can be more easily understood as a function of pres- 
sure with fixed temperature, as in Fig. @] Note that, 
while approaching the BF-DF transition through a pres- 
sure decrease, a is diverging to — oo in the BF side and to 
+oo on the DF side, as discussed in the previous section. 

In the next section we turn our attention to Monte 
Carlo simulations used to study the anomalous behavior 
of the diffusion constant of this model. 
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FIG. 5. Self-Diffusion constant as a function of density at 
fixed temperature. 



IV. MONTE CARLO SIMULATIONS 

Monte Carlo simulations were performed in the canon- 
ical ensemble using jumps to nearest and next-nearest 
neighbour sites. Next-nearest neighbour jumps are nec- 
essary to observe diffusion in one dimensional lattice flu- 
ids since it is impossible for the particles to turn around 
their neighbours (22j . 

Each simulation starts with TV particles placed ran- 
domly along V sites of a linear lattice with periodic 
boundary contourn, at temperature t. The initial con- 
figuration is recorded in a vector X(0) with the position 
of each molecule. On each MC timestep N particles are 
randomly selected to jump and, in the absence of a col- 
lision, the Metropolis algorithm is applied to test the 
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FIG. 6. The location of the maximum (filled circles) and 
minimum (empty circles) values of the self-diffusion constant 
are compared to the line of temperature of maximum density 
(line) in the t vs. p phase diagram. 



acceptance of the movement 23] . After an equilibration 



time T eq the mean square displacements in relation to the 
initial configuration is calculated and recorded at every 
r r timesteps until the simulation time r/ . A number 1Z of 
different simulations is performed to avoid dependences 
on initial conditions. The diffusion constant is then ob- 
tained by interpolating the Einstein equation: 



N 



J2lXk(r) - X k (Q)} 2 ) =2D(t,p)r 



(7) 



\k=l 



■R 



where the average is performed over the 1Z initial condi- 
tions. 

This protocol was used to simulate a system of size 
V = 1000, with the temperature t ranging from 0.30 
to 0.90 and the density in the interval 0.18 < p < 0.9. 
The simulation times used for this system was r/ = 10 7 , 
r eq = 1.2 x 10 6 and r r = 5 x 10 3 ; and a number TZ = 10 
of different initial configurations were performed on each 
point. The error on the diffusion constant was smaller 
than 0.5% in all cases. 



Anomalous Diffusion and Hierarchy of 
anomalies 



The self-diffusion is shown in Fig. [5] as a function of 
density for five different temperatures. At high temper- 
atures, and particularly for t — 0.7 in Fig. [5j the dif- 
fusion decreases with density, as expected in a system 
with simple hard core particles. For temperatures be- 
low t w 0.55 it is possible to observe points of minimum 
(D min ) and maximum (D max ) diffusion, as a function of 
density. Between D m i n and D max the self-diffusion con- 
stant anomalously increase with density, as happens for 



water in a region of the phase diagram [1] . The presence 
of a minimum in diffusion was observed in atomistic Q 
models of liquid water, as well as continuous [IH [24| and 
lattice [25T - I27] models with waterlike behavior but, as far 
as we know, this is the first time that this behavior is 
observed in a one dimensional lattice model. 

In order to locate D m i n and D max points on the neigh- 
bourhood of minimum and maximum diffusion were fit- 
ted with the function: 

D(t, p) = (ap 3 + bp 2 +cp + d) exp (ep 2 +fp + g), (8) 

with a, b, c, d, e, /, and g as numerically adjustable 
parameters. At low temperatures (t < 0.43), two pa- 
rameter sets were used to fit D(t,p), with at least 10 
points around each extrem. For temperatures higher 
than t = 0.43 a single function was used to fit at least 17 
points covering both extrems. Maximum and minimum 
values of diffusion were calculated through the derivative 
of eq. ©. Points of extremum diffusion are compared 
to the exact TMD line in the density vs. temperature 
phase diagram on Fig. [5] Even though not reaching the 
typically lower temperatures of the region anomalous in 
density, the simulations indicate that for less dense states 
(p < 0.6) the line of minimum diffusion occurs at temper- 
atures higher than the TMD line. For higher densities, 
the line of maximum diffusion presents a reentrance and 
it crosses the TMD line, indicating that an strict hierar- 
chy of anomalies is absent on this simplified model (at 
least for e vdw /fhb = 2/3). 

By comparing Fig. [5] and similar results from other lat- 
tice models |25l - l27j one observes an overall tendency of 
overlap between regions of anomalous diffusion and den- 
sity. Nevertheless, the so-called hierarchy of anomalies 
does not always occur in simple hierarchical structure, 
as found in atomistic models of water [||, [2lJ. In two- 
dimensional lattice models with waterlike behavior a hi- 
erarchy of anomalies was found in a reversed order, i.e., 
with the TMD line covering the anomalously diffusive 
region (2(| H3| ■ On the other side, a three dimensional 
lattice model also presents a line of extremum diffusion 
crossing the TMD, as in this work, but with the region of 
diffusion anomaly covering the TMD at low densities [25[ . 
Considering these models and the results presented on 
this paper it is possible to conclude that there is an over- 
all tendency for the concomitant ocurrence of density and 
diffusion anomalies but the existence of an hierarchical 
order between these anomalies depends on dimentional- 
ity and even on possibily on the detailed nature of the 
interactions between molecules. 



V. CONCLUSION 

We introduced a simple one dimentional lattice- 
gas model that reproduces waterlike anomalies on its 
thermodynamics and dynamics through transfer ma- 
trix technique and Monte Carlo simulations. Although 
other models were used to investigate water's peculiar 
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properties[15l ll8H20j . (as far as we know) this is the first 
time this type of model is used to investigate the inter- 
relations between thermodynamic and dynamic waterlike 
anomalies. 

Our results show that the region of thermodynamic 
anomalies is located in the neighbourhood of the transi- 
tion between two structured fluids: a bonded fluid (BF) 
and a dense fluid (DF). The null temperature 'second 
critical point' is linked to this transition and the relation 
between the criticality and the temperature of maximum 
density can be understood from the behavior of the ther- 
mal expansion coefficient in the vicinity of the critical 
point. 

An increase on the diffusion coefficient with density (at 
fixed temperature) was also found in the neighbourhood 
of BF-DF transition. Even though existing an overall 
tendency for the concomitant ocurrence of the regions 
of anomalous diffusion and density, a strict hierarchy of 
anomalies is absent on this model. 

As final remark we mention that the simplicity of the 
current model makes it an ideal prototype for more de- 
tailed investigations on the inter-relations between wa- 
terlike thermodynamic and dynamic anomalies. We 
are now working on analytical approaches to the diffu- 
sion constant of this model using techniques from non- 
equilibrium statistical mechanical. Another approach 
that is currently under investigation involves the calcula- 
tion of the diffusion constant through the Kubo relation 
to understand the connection between memory and dif- 
fusion in more complex topologies (28j . 
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(A3) 



where the sum is made over all lattice states f) = 
{Vo,...,Vl-i} and Q(r)i-i,rn,r) i+l ) = e -W^-i** +1 ). 
Let us now define an auxiliar occupation variable rj[ and 
a new Boltzmann weight as 



,TliTWi- 



(A4) 



ft is not difficult to use P(r]i-i, ?y,, rf t , to write (|A3[) 
as trace of a matrix: 



S(T,L,fi)=H{V L } 



(A5) 



where the elements of matrix V are ordered through bi- 
nary numbers whose bits are occupation variables, and 
are given by: 

[P] m ,^,, v „, = P(r],r]',r]",r]"'). 

In the thermodynamic limit, only the largest eigen- 
value (A) of V contributes to the right hand side of 
eq. (|A5[) . and this equation becomes 
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Appendix A: Transfer matrix technique 

Here we use the transfer matrix technique to calcu- 
late the grand canonical partition function of linear chain 
with Hamiltonian ([lj. More detailed survey of the tech- 
nique can be found on textbooks of statistical mechan- 
ical [13, HH . Since our model does present interactions 
up to three neighbouring lattice sites, we will rewrite the 
Hamiltonian as 



U 



h(r)i-i,rn,r)i + i), 



(Al) 



where periodic boundary conditions are assumed 
implicitly. The three-site interaction Hamiltonian 
h(r]i-i,r)i,r]i + i) is defined as 

h(rj i -i,rii ) r) i+ i) = --e vdw (rji-irji + mVi+i) 

-e hb r)i^i(l - rji)r) i+ i - jUJ^. (A2) 



5 = X L = exp(pPL), (A6) 
where the characteristic polynomial of V is 

A 3 - (1 + za)\ 2 + z{a - b)X + z(b - 1) = 0. (A7) 



By solving together equations (|A6|) and (|A7|) one can 
obtain either P = P(T,n) or /i = /i(T,P). Here we 
follow Bell in ref. [T8(, who noted that z — e 13 ^ can be 
rearranged from (|A7|) as a function of A, a and 6; leaving 
the possibility of calculating the density as 



z /<91n^ 
AL 



dz 
OX 



b 

-i -1 



a , b 



(A-l)[(aA + 6)(A-l) + l] 
(A - f ) 2 (aA + 26) + 3A - 2 



(A8) 



Finally, the same procedure can be used to calculate 
the densities of hydrogen bonds and nearest neighbours 
(per site) as 
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Phb = 



A[3A 2 -2(1 + za)\ + z(a-b)Y 
za 

\{l + bz/\ 2 + z(2\- 1)[A(A- l)]" 2 }' 



(A9) 
(A10) 
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